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I. Introduction 

A djoint solutions of the governing flow equations are becoming increasingly important for the devel- 
opment of efficient analysis and optimization algorithms. A well-known use of the adjoint method is 
gradient-based shape optimization. 1 * " 7 Given an objective function that defines some measure of performance, 
such as the lift and drag functionals, its gradient is computed at a cost that is essentially independent of the 
number of design variables (e.g., geometric parameters that control the shape). Classic aerodynamic applica- 
tions of gradient-based optimization include the design of cruise configurations for transonic and supersonic 
flow, as well as the design of high-lift systems. 

Cartesian-mesh methods 8 " 12 are perhaps the most promising approach for addressing the issues of flow 
solution automation for aerodynamic design problems. In these methods, the discretization of the wetted 
surface is decoupled from that of the volume mesh. This not only enables fast and robust mesh generation 
for geometry of arbitrary complexity, 13, 14 but also facilitates access to geometry modeling and manipulation 
using parametric computer-aided design (CAD). In previous work on Cartesian adjoint solvers, Melvin et 
al 15 developed an adjoint formulation for the TRAN AIR code, 8 which is based on the full-potential equation 
with viscous corrections. More recently, Dadone and Grossman 16, 17 presented an adjoint formulation for the 
two-dimensional Euler equations using a ghost-cell method to enforce the wall boundary conditions. 

In Refs. 18 and 19, we presented an accurate and efficient algorithm for the solution of the adjoint 
Euler equations discretized on Cartesian meshes with embedded, cut-cell boundaries. Novel aspects of the 
algorithm were the computation of surface shape sensitivities for triangulations leased on parametric-C AD 
models and the linearization of the coupling between the surface triangulation and the cut-cells. The accuracy 
of the gradient computation was verified using several three-dimensional test cases, which included design 
variables such as the freestream parameters and the planform shape of an isolated wing. 

The objective of the present work is to extend our, adjoint formulation to problems involving general shape 
changes. Factors under consideration include the computation of mesh sensitivities that provide a reliable 
approximation of the objective function gradient, as well as the computation of surface shape sensitivities 
based on a direct-CAD interface. We present detailed gradient verification studies and then focus on a shape 
optimization problem for an Apollo-like reentry vehicle. The goal of the optimization is to enhance the lift- 
to-drag ratio of the capsule by modifying the shape of its heat-shield in conjunction with a center-of-gravity 
(c.g.) offset. This multipoint and multi-objective optimization problem is used to demonstrate the overall 
effectivness of the Cartesian adjoint method for addressing the issues of complex aerodynamic design. 
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II, Optimization Problem 


The aerodynamic optimization problem we consider in this work consists of determining values of design 
variables that minimize a given objective function 

min J{X,Q) (1) 

where J represents a scalar objective function defined by a surface integral, for example lift or drag, X de- 
notes a scalar design variable, for example a shape parameter of the wetted surface, and Q — [p, pu, pv, pw , pE) T 
denotes the continuous flow variables. The flow variables are forced to satisfy the steady-state three- 
dimensional Euler equations within a feasible region of the design space O 

X(X,Q) = 0 VleO (2) 

which implicitly defines Q = f (X). Our goal is to compute a reliable approximation to the objective function 
gradient dj jdX . We proceed by using the discrete adjoint method, where the governing equations, Eqs. 1 
and 2, are first discretized and then differentiated. 

A. Discrete Flow Equations and Solution Method 

The discretization of the Euler equations uses a second-order accurate finite-volume method on a multilevel 
Cartesian mesh with embedded boundaries. The . mesh consists of regular hexahedral cells, except for a 
layer of body-intersecting cells, or cut-cells, that are arbitrary polyhedra adjacent to the boundaries. A 
cell-centered approach is used, where the control volumes correspond to the mesh cells and the cell-averaged . 
value of Q, denoted by Q } is located at the centroid of each cell. The spatial discretization uses the flux-vector 
splitting approach of van Leer. 20 The boundary conditions are enforced weakly by appropriate modifications 
of the reconstructed state and boundary flux. The resulting discrete system of equations is given by 

R(Q 1 M i X) = Q (3) 

where Q ~ [Qi, Q 2> * • * , Qn] t is the discrete solution vector for all N cells of a given mesh M, and R is the 
flux residual vector. Steady-state flow solutions are obtained using a five-stage Runge-Kutta scheme with 
local time stepping, multigrid, and a highly-scalable domain decomposition scheme for parallel computing. 
For further details on the spatial discretization and flow solution, see Aftosmis et a/. 14,21,22 and Berger ei 
aZ . 23, 24 

Note that design variables appear directly in Eq. 3 only when they involve parameters that do not change 
the computational domain, such as the Mach number, angle of incidence, and side-slip angle. The influence 
of shape design variables on the residuals is implicit via the computational mesh 

M == f[T(X)) " (4) 

where T denotes a triangulation for the surface model. The functional dependence of the triangulation on 
the design variables is determined by the geometry parameterization scheme and is explained further in 
Sec. III. 


B. Discrete Adjoint Method 

Combining Eqs. 3 and 4 and differentiating about a steady-state solution Q gives the gradient of the discrete 
objective function J(X , M, Q) 


d dJdM dT 

AX ~ dX + dM&T dX 


l 






8R , d_R d_M_ df 
dX 1 m df dX 


A 


7 


where the vector ■>!> represents adjoint variables given by the adjoint equation 


8R' 


dj ■ 


dQ V> ~~d§ 


( 5 ) 


( 6 ) 
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We emphasize that the adjoint equation is independent of the design variables, i.e., the partial derivatives 
with respect to the flow variables are evaluated at constant X. Hence, for problems with many design 
variables and one objective function, the use of the adjoint method is attractive since only one solution of 
the large linear system of Eq. 6 is required. 

The solution algorithm for the adjoint equation leverages the Runge-Kutta time-marching scheme and 
the parallel multigrid method of the flow solver. The algorithm is implemented using the duality-preserving 
approach, 25 such that the asymptotic convergence rate of the discrete adjoint is identical to that of the 
flow solver. The matrix- vector products associated with the flow- Jacobian matrix, left side of Eq. 6, are 
computed on-the-fly using a two-pass strategy over the faces of the mesh. The flow- Jacobian matrix, as well 
as the term dJ/dQ in Eq. 6, are derived by hand, where we neglect the linearization of the limiter function 
used in the solution reconstruction procedure. Overall, the CPU time per iteration and memory usage of the 
adjoint solver are roughly equivalent to the flow solver. For further details on the adjoint formulation and 
solution procedure see Ref. 18. The computation of partial derivative terms dJ/dX and dR/dX in Eq. 5 
is straightforward, since these terms do not involve derivatives of the surface shape. The remaining partial 
derivative terms in Eq. 5, labeled as A and B, represent the differentiation of the objective function and 
residual equations with respect to design variables that alter the surface shape. An accurate computation of 
these terms is presented in Ref. 19 and an overview is given in the following section. 

III. Objective and Residual Shape Sensitivities 

The differentiation of the objective function and residual equations with respect to shape design vari- 
ables involves the interaction of the surface deformations with the triangulation and the volume mesh. An 
advantage of our Cartesian approach is that the computation of the surface shape sensitivities, i.e., the term 
dT'/dX in Eq. 5, is cleanly isolated from the residual equations. This is because there is no prescribed con- 
nectivity between the volume mesh and the triangulation, which is particularly well-suited for CAD-based 
geometry control. In subsection A, we outline our approach to geometry manipulation and surface trian- 
gulation using a direct CAD interface, and then explain the computation of surface shape sensitivities. In 
subsection B, we address the binding of the surface to the volume mesh to obtain sensitivities in the cut-cell 
layer. 

A. Shape Sensitivity of Surface Triangulation 

Our primary method for surface modeling, control, and triangulation is based on the Computational Analysis 
and PRogramming Interface (CAPRI) developed by Haimes et al 26 ’ 27 CAPRI provides a unified and 
direct access to most parametric CAD systems. The ability to control CAD solid models is accomplished 
by exposing the master-model feature tree of parametric parts and assemblies. Design variables can be 
associated directly with the exposed parameters, or indirectly with splines, such as airfoil sections that are 
lofted to define a wing. 

Upon a successful regeneration of the CAD-model, CAPRI provides an associated water-tight triangu- 
lation for each component. The Cartesian position of each vertex in the triangulation, for each face of the 
CAD model, is given by its parametric values 

(x,y,z) = f(u,v) (7) 

Since analytic differentiation of the geometry constructors within the CAD system is presently not possible, 
we rely on a centered- difference approximation to determine the derivative of each vertex of the triangu- 
lation with respect to the design variables. The first step of this procedure involves the regeneration and 
triangulation of the CAD model to reflect the current values of the design variables, thereby obtaining a 
baseline model. For each CAD face, we normalize the (u, v) vertex values such that v 6 [0, 1]. The next 
step involves two additional model regenerations and triangulations for each design variable, which corre- 
spond to the plus and minus perturbations. To form the finite-difference approximation, we interrogate the 
perturbed models using the normalized (u,v) values re-scaled by their new (u, v) range. Hence, we obtain a 
one-to-one correspondence between vertex coordinates in the physical space for the baseline model and its 
perturbations. A tacit assumption in this procedure is that the face topology of the baseline and perturbed 
geometries remains the same. We check this condition after each regeneration of the model and attempt to 
reduce the stepsize or use one-sided differences if necessary. 
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The finite-difference procedure is efficient since we query the triangulation in the parameter space, i.e., 
we only require forward solutions of Eq. 7. An alternate approach is the use of “snap” operations in the 
physical space to locate nearest vertex coordinates on the baseline and perturbed surfaces. This procedure 
is less efficient since it involves an inverse solution of Eq. 7. Furthermore, due to internal tolerances of the 
CAD system, snap operations may introduce excessive noise into the evaluation of the derivative. 


B. Sensitivity of Cut-Cells to Shape Deformations 

The cut-cell generation procedure involves computing the intersection of the surface triangulation with the 
faces of Cartesian hexahedr a. Details of this procedure are described in Ref. 13. The sensitivity of a cut-cell 
to a shape deformation is obtained by a linearization of a geometry constructor that closely approximates the 
mesh generation process. We approximate the cut-cell geometry using the intersection of surface triangles 
with only the Cartesian edges and neglect contributions from the faces. This is illustrated in Fig. 1, where a 
Cartesian edge de of a hexahedron intersects triangle (a,b, c), creating a pierce point P. We use bold type 
to represent Cartesian position vectors. 

The linearization of the pierce point location with respect to 
a shape perturbation could be obtained directly from the CAD 
system using the approach described in the previous section. 

This would involve the use of “snap” operations to locate the 
pierce point in the (u : v) space of the appropriate CAD face. 

Alternatively, the pierce point location can be defined using a 
set of parametric equations for the triangle (a, b, c) and the 
line segment de, which is the approach we pursue in this work. 

The location of P along the segment de is given by 


P — d.+ s* • D 


( 8 ) 



where D ~ e — d and s* = /(a, b, c, d, e) represents the para- 
metric value of the intersection point. The expression for s* 
is derived in Ref. 13. The linearization of the pierce point is 
given by 

9P _ 

dX ~ dX D 

Note that only one component of D is non-zero, since the Cartesian mesh points are coordinate aligned. 
Furthermore, the points d and e are independent of X, i.e., the Cartesian mesh is rigid, and consequently 
the linearization of s* is given by 


Figure 1. Intersection of a triangle with a 
Cartesian edge defining a pierce point P 


( 9 ) 


_ ds* da. ds* db ds* dc 
dX = da dX + 'dbdX + dc dX 


( 10 ) 


This linearization represents the change in the location of P along the edge de exclusively as a function of 
shape sensitivities at the vertices (a,b,c). The decoupling of the surface sensitivities from the volume mesh 
provides a straightforward interface to any geometry parametrization tool. 

The linearization of the pierce points is used to determine all the required geometric sensitivities in the 
linearization of the residual equations. First-order spatial discretization requires the derivative of the face 
areas and normals, while for second-order discretization, we require additional linearizations for the location 
of face centroids, volume centroids, and the reconstruction distances used in the evaluation of flow-solution 
gradients. As a result, a perturbation of the surface shape influences the residuals not only in the cut-cells, 
but also in the first and second nearest neighbouring cells. 

The computation of the surface normal is based on a planar approximation to the variation of the trian- 
gulation wdthin each cut-cell. This agglomerated normal vector is computed using the divergence theorem, 
which requires the geometric closure of each cut-cell 

D = D A i • = 0 ( n ) 

jeVi jev, 

where n denotes an area scaled normal vector. The sum is performed over the Cartesian faces of the cell, 
which must equal the agglomerated normal vector. Since geometric closure of a cut-cell should be satisfied 


4 of 10 


American Institute of Aeronautics and Astronautics 



for any shape perturbation, we linearize Eq. 11 to obtain 

. ( 12 ) 

: j€Vi ' 

The derivative of the Cartesian face areas is computed based on the linearization of the pierce points. The 
pierce points, in conjunction with the corner points of the cut-cell, form polygons on each Cartesian face. The 
areas of the polygons are computed by subdividing each polygon into triangles. Central to this procedure is 
the linearization of an area of a triangle given by the cross-product of two vectors. 

The final paper will include a discussion of the computation of volume mesh sensitivities used to obtain 
accurate gradients. 

TV. Results and Discussion 

We present several three-dimensional test cases to verify the accuracy of the cut-cell linearization and 
gradient computations. The Pro /EN GINEER® Wildfire CAD system is used to create all geometry models. 

A. NACA 0012 Wing 

This test case is used to verify gradient accuracy and involves an isolated wing at transonic flow conditions, 
Moo = 0,84 and a = 3.06°. The wing geometry is constructed from a generic CAD model with linearly lofted 
NACA 0012 root and tip sections, a taper ratio of 0.7, aspect ratio of 11, sweep of 15°, and root and tip-twist 
angles of 3° and —1.5°, respectively. The volume mesh contains 499,716 cells and uses a symmetry plane. 
Convergence to steady-state is achieved using 64 processors 8 -, a 4-level W-cycle multigrid with one pre- and 
one post-smoothing pass, and a CFL number of 1.2. Partial updates of the flow gradients are also used, i.e., 
the flow gradients are updated only on the first stage of the RK5 scheme. The second-order accurate spatial 
discretization uses the van Leer limiter. The Mach number contours of the flow solution are shown in Fig. 2. 
The objective function represents a lift-constrained drag minimization problem 

J -( 1 -|) 2+0 ' 01 ( 1 ’i) 2 (13) 

where == 0.327 denotes the initial lift coefficient and — 0.01. The initial drag coefficient is 0.0487. 
The design variables are the angle of incidence and the tip-twist angle of the wing. 

Convergence of the flow and adjoint equations are shown in Fig. 3. Approximately the first 50 multigrid 
cycles of the flow solution and the first 75 cycles of the adjoint solution correspond to full-multigrid startup. 
As shown in Fig. 3(b), the startup time is almost negligible. Furthermore, the convergence characteristics of 
both solvers are quite similar, with a reduction in residual of over six orders of magnitude in 60 seconds. Note 
that the flow solver takes advantage of specific compiler directives for code optimization. These directives 
are not yet implemented in the adjoint code. 

Objective function gradient accuracy is presented in Table 1 for both design variables. The agreement 
between the adjoint and centered-difference gradient values is good. The small differences are attributed to 
the constant limiter assumption in the linearization of the discrete flow equations, as well as the dependence 
of the finite-difference gradient on the stepsize due to non-smooth changes in the shock location. The 
interpretation of the descent direction, which is the negative of the gradient vector, reflects the fact that the 
objective function is dominated by the drag term. Hence, a reduction in the angle of attack and an increase 
in the washout of the wing should lower the shock strength and therefore reduce the drag penalty. 

The final paper will include additional design variables , as well as an optimization example for this 
problem . 

B. Reentry Capsule 

The second design example targets the optimization of a heat-shield shape for a reentry capsule. The 
objective of the optimization is to enhance the L/D of the capsule, thereby improving its trajectory control 
and landing-site selection, as well as reducing the reentry load factor and heat rates. The nominal reentry 

a Intel® 1.5 GHz IA-64 Itanium 2 processor 
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Multigrid Cycles WallciockTime (s) 

(a) Multigrid Cycles (4- level multigrid with full- (b) Wallclock Time (64 Intel® 1.5 GHz Itanium 2 
multigrid startup) processors, mesh size 499,716 cells) 

Figure 3. Convergence histories for the NACA 0012 Wing (Moo — 0.84, a — 3.06°) 


Table 1. Gradient accuracy for the 
NACA 0012 Wing 


Design 

Variable 

Adjoint 

Finite-Difference 

.a 

0.172276 

0.177888 

Tip Twist 

0.062750 

0.063810 
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trajectory assumes a constant LjD, winch leads to a multipoint optimization problem. The pitch stability 
of the capsule, i.e., trim and negative C ma slope, is ensured by introducing additional penalty terms in the 
objective function. 

The baseline configuration is shown in Fig. 4(a) and is derived from the Apollo capsule shape . 28 The 
heat-shield is modeled as a two-directional blended surface. Its shape is controlled by the center-line (red 
line in Fig. 4(a)) and surface tangency conditions at the shoulder. The center-line is parameterized using 
a cubic B-spline, as shown in Fig. 4(b). The design variables are associated with three control points 
indicated in Fig. 4(b), each having a single degree-of-freedom in the vertical direction. Additional structural 
constraints, such as heat-shield length, convex shape, and radius of curvature, can be introduced to account 
for aerothermal factors. • 



(a) Solid model and coordinate system 


— — Heat Shield 
□ Control Points 

b Design Variables 



(b) Shape parameterization of heat shield center-line (normalized by cap- 
sule diameter d — 5.5m) 


Figure 4. Parametric- CAD model of capsule shape 

We consider a two-point optimization problem where the design Mach numbers are 10 and 25. The 
corresponding values of 7 , the ratio of specific heats, are 1.231 and 1.125, respectively. The target value of 
LjD at each design point is set to 0.4. This value is based on the aerodynamic characteristics of the Apollo 
capsule, which attained an LjD of 0.3 using a c.g. offset. We use a similar c.g. offset in the present study, 
given by x/d = 0.124, y/d = 0, and zjd— 0.03698, where d represents the capsule diameter set to 5.5m. 
The objective function is given by 

j = 0.1 (| -0.4Y 2 + Cl + O.lD-0 A 
- ^ ^ J M— 10 L • 
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The angle of incidence at each design point is used to enforce the pitching-moment constraint, resulting in 
a total of five design variables. The gradient is computed using finite-differences. 

Figures 5 and 6 shown the baseline and optimal designs at the two operating conditions. The shape 
modifications required to achieve the target performance are" relatively minor. The ^objective function is 
satisfied within 15 design iterations and the gradient is reduced* by roughly 2.5 orders of magnitude, as 
shown in Pig. 7, Figure 8 shows that value of L/D improved not only for the design Mach numbers, but 
also at off-design points. The variations in L/D for the final design, as shown in Fig. 8, could be reduced by 
introducing an additional operating condition in the optimization problem. 

In the final paper, the adjoint method will be used to compute the gradient We will also discuss the pitch 
slope constraint and investigate c.g . sensitivities . Additional aerodynamic coefficients will be presented and 
discussed. 



Figure 5. Pressure contours at M = 10 (Baseline atrim. — 156.1 deg., Optimal Shape a-trim = 155.7 deg.) 
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